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Abstract. We study the complexity of computing the real solutions of 
a bivariate polynomial system using the recently proposed algorithm Bl- 
SOLVE [3]. Bisolve is a classical elimination method which first projects 
the solutions of a system onto the x- and j/-axes and, then, selects the 
actual solutions from the so induced candidate set. However, unlike sim- 
ilar algorithms, Bisolve requires no genericity assumption on the input 
nor it needs any change of the coordinate system. Furthermore, exten- 
sive benchmarks from [3] confirm that the algorithm outperforms state 
of the art approaches by a large factor. In this work, we show that, for 
two polynomials /, g £ Z[x, y] of total degree at most n with integer co- 
efficients bounded by 2 T , Bisolve computes isolating boxes for all real 
solutions of the system / = g = using 0(n a r 2 ) bit operations^, thereby 
improving the previous record bound by a factor of at least n 2 . 



1 Introduction 

Systems of polynomial equations naturally arise in many fields of science and 
engineering. In computational geometry and computer graphics, there is a par- 
ticular interest in the study of polynomial systems in two or three variables. For 
instance, all existing exact and complete algorithms for computing the topology 
or arrangements of algebraic curves |6ll3j (and surfaces [3]) are crucially based 
on determining the critical points, which are in turn the solutions of a bivariate 
polynomial system. In this work, we investigate in the bit complexity analysis 
of the recently proposed algorithm Bisolve [3] to isolate the real solutions of a 
bivariate polynomial system 

f(x,y)= /<r r V °> 9{x,y)= J2 £^y=0, (1.1) 

i,j£N:i+j<n i,j£N:i+j<n 

where f,g£ Z[x,y] are polynomials of total degree at most n and with integer 
coefficients bounded by 2 r , r £ N. For short, we will also write that / and g 
have magnitude (n, r). Henceforth, we assume that / and g share no common 
non-trivial factor in Z[x, y]\Z which, due to Bezout's Theorem, is equivalent to 
the existence of finitely many (complex) solutions of Bisolve computes a 

set of disjoint boxes Bk C K 2 such that the union of all Bk contains 

Vm := {(x,y) € R 2 \f(x,y) = g{x,y) = 0}, 



1 O indicates that polylogarithmic factors in r and n are omitted. 



the set of all real solutions of and each Bk is isolating, that is, it contains 

exactly one solution. We show that to achieve the latter task, Bisolve demands 
for 0(n S T 2 ) bit operations, improving the previous record bound by a factor of 
n 2 or more (depending on whether t or n is dominating) . Our analysis uses two 
recent results on isolating [17] and refining [14] the real roots of a univariate 
polynomial. Yet, we remark that the so obtained complexity bound for Bisolve 
is not merely due to the improved complexity bounds for the root isolation and 
refinement step, but also due to the effectiveness of the novel inclusion predicate 
used in the validation step of Bisolve to certify or discard candidate solutions. 
Wc remark that, in all previous algorithms to compute Vr, the computation of 
the common roots of f(a,y) and g(a,y) was achieved by considering a corre- 
sponding signed remainder sequence, with a the projection of a solution in R 2 
onto the real axis. In practice, the latter computation often turns out to be a 
bottleneck of the overall approach. Bisolve does not follow this approach which 
has proven to be favorable in practice. Our analysis also shows that the final 
validation step is less hard than isolating the roots of the resultant polynomial, 
which seems to be the bottleneck. In addition, wc would like to emphasize that 
BiSOLVE is not a "galactic" algorithm just designed to achieve good complexity 
bounds. In contrast, the high performance of Bisolve in practice is confirmed by 
many experiments; we refer to [3] for an extensive comparison of the algorithm 
with other state of the art approaches such as Lgp [7] or Maple's Isolate. 

Related work. An early result on the complexity analysis appears in [12j . where 
the closely related problem of computing the topology of an algebraic curve 
is considered. The authors analyze the algorithm Top and derive a complexity 
bound of (D(N 1A ) bit operations, with N = max(n, r). Another work [5] presents 
several methods to solve a bivariate polynomial system. The first method, Grid, 
first projects the solutions onto orthogonal axes and, then, matches them by 
means of a Sign_AT procedure which uses evaluation of a signed remainder se- 
quence of / and g. The complexity of Grid is bounded by C(iV 14 ) bit operations, 
where the overall cost is dominated by the SiGN_AT operations. The second ap- 
proach, M_RUR, is based on subresultants and RUR (rational univariate repre- 
sentation) techniques and achieves a bit complexity O(n 10 (n 2 + r 2 )) = 0(N 12 ). 
The third approach, G_RUR, achieves the same complexity as M_RUR but re- 
lies on the computation of the GCD of the square- free parts of /(a, y) and g(a, y). 
It should not be revealed that, using the improved complexity bounds for the 
univariate root isolation, one may also obtain improvements of the correspond- 
ing complexity bounds for M_RUR (only in a sheared system) and G_RUR. 
However, it seems that the so obtained bounds are still weaker than the bound 
achieved by Bisolve. Namely, the dependence on n in the final steps of M_RUR 
and G_RUR is considerably larger (i.e., by a factor n 2 when n is dominating). 
We further remark that the complexity analysis of G_RUR is based on the 
modular algorithm in |19j whose complexity is only given with respect to aver- 
age running time (cf. Section 3.2 in [T5]). In his dissertation, M. Kerber describes 
randomized algorithms to analyze the topology of a single algebraic curve and 
to compute arrangements of such curves. A detailed analysis of the "curve-pair 



analysis" which solves the subproblem of finding the solutions of a bivariatc 
system shows that the corresponding bit complexity is bounded by an expected 
number of O(n 10 (n + r) 2 ) bit operations; see [T31 Section 3.3.4]. 

Outline. Section [2] introduces some notations which are used throughout the ar- 
gument. In Section [31 we briefly review the algorithm Bisolve omitting some 
technical details and filtering techniques to keep the presentation simple. The 
complexity analysis is given in Section [4j We analyze the three main steps of 
the algorithm separately and, then, combine the results yielding the overall com- 
plexity. Finally, in Section [SJ we give some concluding remarks. 

2 Setting 

We express the input polynomials / and g in as univariate polynomials in 
x and y of degrees n x and n y , respectively: 

f(x,y) = £jf \y)x i = g/^W, 9{x,y) = £ ff^ (l/)^ = X> Cw W. 

i=0 i=0 i=0 i=0 

where Z^, c/^ G Z[ x] and ft , gf g Z[y]. Throughout the paper, it is assumed 
that n x , n y < n. We denote the Sylvester's matrix associated with / and <? by 
S (y) = sM(f,g) whose entries are the coefficients {f} y) } and see [IS 

p. 286] for the definition. The resultant R^ = tes(f,g,y) g Z[x] of / and 3 
with respect to the variable y is the determinant of S^ v '. By analogy, = 
res(/, <?,x) defines the resultant with respect to the variable x, and S^ x '(f,g) is 
the associated Sylvester's matrix with entries {f- } and {g^}. If this causes 
no ambiguity, we also write R omitting the variable index and by R* the square- 
free part of R. For an interval I = (a, b) C R, wj := b — a denotes the width, 
mj := (a + b)/2 the center and rj := (b — a)/2 the radius of I. A disc in C is 
denoted by A = A r (m), where m g C defines the center of A and r G R + its 
radius. For a (not necessarily square-free) polynomial F(x) = X]™=o ^ iX% e 
with distinct roots zi . . . z m g C, the separation <Ti := a(z.i,F) of a root 2j is 
defined as the minimal distance of to any root Zj 7^ z*. The root separation 
cr(F) of i* 1 is defined as the minimum of all <7j, and -E'(.F') = X)2=i l°gcr(z,;, F) -1 . 
We finally denote r(F) := max^ |jZj| the maximal absolute value of all Zi. 

3 Review of the algorithm 

In this section, we recall the algorithm Bisolve to make the paper self-contained; 
for further details and filtering techniques used in the actual realization, we re- 
fer to [3]. At the highest level, Bisolve comprises three subroutines which we 
consider separately. 



Project: We begin with projecting the complex solutions of onto the x- 
and y-axes. In other words, we consider the two sets: 

V^ x) := {x G C\3y G C A f(x, y) = g(x, y) = 0} 
V™ := {y G C|3x G C A f(x, y) = g(x, y) = 0} 

and compute their restrictions :— fll and := DR to the real 
values. The real solutions Vr of are then contained in the product 

C := x V^' y) C R 2 , (3.1) 

which we denote the set of candidate solutions for . To obtain the projection 
sets Vj x) and V^, we first compute the resultants R^ and RS X \ respectively, 
and extract the square-free part R* of either polynomial (R — R^ or R — R^ 
for short). Next, we isolate the real roots a, of R* using RIsolate, a recently 
proposed Descartes method [17] mentioned in the introduction. RIsolate is an 
exact method which works with approximations of the coefficients of a polyno- 
mial and achieves a significantly better bit complexity compared to the classical 
Descartes algorithm [3J Remark 10.52]. Moreover, in contrast to asymptotically 
fast algorithms such as those proposed by Pan or Schonhage (see [TB] for an 
overview), RIsolate is practical and easy to implement; see [T5] for an imple- 
mentation within Mathematica. The isolating intervals I := 1(a) returned 
by Msolate satisfy 

cr(a,i?)/(4degi?*) 2 < wj < 2(dcgiT) ■ a(a,R), (3.2) 

which will be important for us in the context of our complexity analysis. This 
concludes the first step of the algorithm. 

Separate: In this step, the real roots of R are further separated from the com- 
plex ones. For each root a, we refine a corresponding isolating interval / := 1(a) 
such that the disc ^^(m/) does not contain any root of R except a. RIsolate 
returns an interval I which fulfills the inequality in (|3.2[) . Thus, after performing 
[log(16 dcg(i?*))] bisection steps to refine /, the disc A% ri (mj) isolates a as well. 
Then, wc compute 

LB(a) := 2- 2dc ^ R \R(m I - 2n)\, (3.3) 

which constitutes a lower bound for |i?(x)| on the boundary of A(a) := Ai ri (mi), 
that is, \R(x)\ > LB (a) for all x G dA(a); see [3 Thm. 3.2] for a proof. We 
remark that in the initial description of the algorithm (as well as in the real- 
ization), we consider a different predicate to ensure that Ag ri (mi) is isolating. 
The latter predicate is based on a local Taylor expansion of R* at mi, docs 
not require (|3.2[) and shows a slightly more adaptive behavior in practice. It is 
rather straight forward to prove that also using this method the overall approach 
achieves the claimed bit complexity bound. However, for the sake of simplicity, 
we decided to exploit the inequality in p.2[) for our complexity analysis. 



To sum up, at the end of Separate, for each real root a of R^ (and /3 
of R&), we have an isolating interval 1(a) (and 1(0)) and an isolating disc 
A(a) := Z\ 2r/(ci) (m 7 ( Q )) (and A(0)). Then, each real solution of system is 
contained in a polydisc A(a, 0) :— A(a) x A(0) C C 2 , and each of these polydiscs 
contains at most one solution. In addition, for each point (a;, y) on the boundary 
of a polydisc A(a,0), wc have \R {y) (x)\ > LB (a) or \R {x) (y)\ > LB(j3). 

Validate: The goal of this final stage is to determine all candidates (a, 0) 
which are actually solutions of (jl.ip and to exclude the remaining ones. Again, 
to facilitate the complexity analysis, we assume that the actual solutions are 
chosen exclusively based on the inclusion test outlined below. We remark that 
the efficiency of the actual implementation is further due to a series of filtering 
techniques to rapidly exclude the majority of candidates such as tests based on 
interval arithmetic or a bitstream Descartes algorithm |10) . 

In Separate, we have already computed lower bounds LB(a) and LB(0) for 
the values of \R^ \ and \R^ \ at the boundaries of A(a) and A(0), respectively. 
We now rewrite R^ in terms of cofactors vS v ' and (see pT] p. 287]): 

R (y) {x) = u (y) (aJ) y )f( X) y) + v (y) {x> y)g ^ y) (3.4) 

where and are determinants of "Sylvester-like" matrices U^ y > and V^ v \ 
These matrices are obtained from S^ v \f,g) by replacing the last column with 
vectors (y n y^ 1 . . . y 1 . . . 0) T and (0 ... . . . y 1) T of appropriate size, 

respectively. Next, we compute upper bounds UB(a, 0u^ v ') and UB(a, v^) 
for \u( yS> \ and \v^\ on A(a,0), respectively. These bounds can be obtained by 
bounding the absolute values of the entries in U^- y > and V^- y ' and then applying 
Hadamard's inequality to U^ v > and V^ v \ We remark that the computation of 
the latter bounds is done without explicitly computing the cofactors (which are 
typically very large expressions). Cofactor polynomials u^ x \ v^ x ' and respective 
upper bounds UB(a,0uW), UB(a,0v^) are defined in an analogous way for 
the resultant polynomial R^ x > . The inclusion test based on a homotopy argument 
is now formulated as follows (see [3j Thm. 4] for a proof): 

Theorem 1. If there exists an (xo,yo) S A(a, 0) with 

UB(a,0u^) ■ \f(x o ,yo)\ + UB(a,0v^) ■ \g(x ,y )\ < LB (a), (3.5) 

UB(a,0u^) ■ \f(x ,y )\ + UB(a,0v^) ■ \g(x ,y )\ < LB(0), (3.6) 

then A(a, 0) contains a solution of il.l]) and, thus, f(a,0) = 0. 

All candidate solutions (a, 0) £ C are now treated as follows: Let B(a,0) = 
1(a) x 1(0) C R 2 be the corresponding candidate box. Each candidate box is then 
refined with the quadratic interval refinement (QIR for short) method from [14] 
until either we can ensure that /(a, 0) ^ or g(a, 0)^=0 based on interval arith- 
metic on B(a,0) or, for an arbitrary point (xo,yo) & B(a,0), the inequalities 
(|3.5I) and (|3.6I) arc fulfilled. In the latter case, Theorem [T] guarantees that (a, 0) 
is a solution of We refer to Section T4.3I for the details of the evaluation 

using interval arithmetic. 



4 Complexity analysis 



Throughout the analysis, we assume that the multiplication of two integers is 
always done in asymptotically fast way. In other words, the bit complexity to 
multiply two k-bit integers is assumed to be M(k) = 0(k log k log log k) = O(k). 

4.1 Project 

Computing the resultant polynomials R = R^ (and R = R^) using the subre- 
sultant PRS algorithm needs (D(n 7 (log n + r)) bit operations. R has magnitude 

(n 2 ,0(n(logn + T))), (4.1) 

see P21 Trims. 2.4.14-17], where R = Sreso(/, g,y), the 0-th subresultant of / 
and g with respect to y. Next, we compute the square-free part R* of R. This 
can be done with 0(n 5 (r + logn)) bit operations (see [13j Thm. 2.4.21]), and 
R* is of magnitude 

(n 2 ,0(n(n + T))). (4.2) 

Finally, the real roots of R* arc isolated using Msolate as outlined in Section[3j 
The complexity of the root isolation is summarized below (see [TTJ Thm. 18]): 

Theorem 2. Let F be a square-free polynomial of magnitude (d,fj,). Then, iso- 
lating the real roots of F demands for 0(d(S(F) + dlog r(F)) 2 ) bit operations, 
with £(F) and r(F) as defined in Section^ 

We aim to apply TheoremHto F := R*. Obviously, r(R*) = r(R) and r(R*) = 
r(R), and a bound on S(R) is provided by the following theorem: 

Theorem 3. For a (not necessarily square-free) polynomial F(x) with integer 
coefficients of magnitude (d, /i) and distinct complex roots zi,...,z r (r < d), it 
holds that E(F) = 0(d/ilog(d/x)). 

Proof. The proof essentially follows the same lines as the proof in [T71 Ap- 
pendix 6.2], with the exception that we use the Davenport-Mahler bound for 
non square-free polynomials. We refer to Appendix [X] for a full argument. 

Now, by Theorem [3] and P~Tj) . it follows that: 

S(R*) = 0(n 2 ■ n{r + log n) • log(n 3 r)) = d(n 3 r). (4.3) 

Observe that, due to Cauchy's bound [2 Corollary 10.4], the absolute values of 
the roots of R (and R*) are bounded by a 

q _ 20(n(logn+r))^ (^A) 

Hence, from (|4.3p . (|4.4[) and Theorem [5J we conclude that the complexity of 
isolating the real roots of R* is bounded by: 

6(n 2 (£(R*) +n 2 log B) 2 ) = 6(n 2 (n 3 T + n 3 r) 2 ) = d(n 8 r 2 ), (4.5) 

which determines the complexity of Project. 



4.2 Separate 



Let a be a real root of one of the resultant polynomials R and / = 1(a) a 
corresponding isolating interval returned by RIsolate. In our description of 
Separate, we have already seen that I = 1(a) is refined by means of O(logn) 
bisection steps to guarantee that the disc A$ ri (mi) isolates a. The bit complexity 
of one bisection step for a polynomial F of magnitude (d, /i) and an interval 
of bitsize 9 is (D(d(fj, + d9)) because this is essentially evaluation of F at the 
endpoints of /; see [131 Prop. 2.5.1] @ The maximal bitsize of the interval / = 1(a) 
is essentially determined by the separation of a. Namely, / is refined to a width 
larger than a(a, i?)/(1024n 6 ) because we started with an interval of width larger 
than cr(a, i?)/(16n 4 ) (see (|3.2j) for a bound on the size of the initial isolating 
interval) and then performed less than 32n 2 further bisection steps. Hence, we 
obtain 6(a) — 0(\ogB + \og(l/a(a,R)) + logn) for the bitsize of 1(a), where 
B is the root bound for R, see (|4.4[) . It follows that the complexity to refine the 
isolating intervals for the real roots z\, . . . , zj~ of R* in Separate is bounded by 
((d, /j,) = (n 2 , (D(n(n + t))) a bound on the magnitude of R*): 

k k 

J2O(d(n + d0( Zl ))) =J20(n\ - n i loga(z l ,R)) = 6(n 7 r), (4.6) 

i=l i=l 

where we use that - YT i= i logcrOi, R) < S(R*) + n 2 log(2B) = 6(n 3 r) because 
each root Zj of R* not which is not considered in the left sum has separation less 
than 2B = d(nr) according to (O and S(R*) = 6(n 3 r) according to (|43)l . 

4.3 Validate 

Estimating lower and upper bounds In the final stage, we have a set of 
candidate solutions C and corresponding disjoint polydiscs A(a, (3) := A(a) x 
Z\(/3) C C 2 . Each of the polydiscs contains at most one solution of (jl.ip (namely, 
(a, /?)). It remains to determine the actual solutions based on the inclusion test 
from Theorem [T] and to exclude the other candidates by means of interval arith- 
metic. We split the complexity analysis in two parts. In the first part, we derive 
a lower bound for the value LB(a) as well as an upper bound for the values 
UB(a, f),iS v >) and UB(a, (3, v^ ) needed by the inclusion predicate. In the sec- 
ond part, we estimate how good each candidate (a, /3) needs to be approximated 
in order to certify or it as a solution or to discard it. 

Estimating the lower bounds. In this section, we derive a lower bound for LB (a) 
which is in turn a lower bound for \R^ \ on the boundary of A(a); sec (|3.3|) for 
the definition of LB(a). By analogy, we obtain a similar bound also for LB(fi) 7 
the lower bound for \R( X > \ on the boundary of A(j3). Remark that the polynomial 
R = R( y *> might have multiple roots. The key idea to keep the bounds tight is to 
exploit the fact that the separation of a root actually depends on its multiplicity. 
We begin with the following auxiliary lemmas. 

2 The bitsize of an interval is defined as maximal bitsize of its rational endpoints. 



Lemma 1. J7J Thru. AB and Thru. 1] Let F be an integer polynomial of degree 
d > 2. Suppose that a is a zero of F{x) of order s and j3 a zero of F{x) of order 
t. If a ^ P and t> s, then 

\a-0\> 2- d/t {d + iyd/t-3d/(2st)^-2d/(st) max ^ max { l5 j^}. 

Lemma 2. Lem. A. 7]. Let F be a non-constant integer polynomial of degree 
d. Let £ be a complex number and a be the root of F(x) which is closest to £ (i.e. 
|£ — a\ is minimal for all roots a of F). Then, denoting by s the multiplicity of 
a as root of F(x), we have: 

\i-a\ s < d d+3d/(2s) \F\ 2 W s -V\F(Z)\. 

Suppose a — zi is a root with multiplicity m, of R and let / = I(zi) be the 
corresponding isolating interval for Zi obtained in Separate. In Section 14.21 
we have already shown that rj > cr(a, i?)/(1024n 6 ) and, thus, the distance of 
z := mi — 2rj to is larger than a (a, R)/(1024n 6 ). Let zj 7^ Zi be another 
root of R closest to Zi (i.e., \zt — Zj\ = o~(zi,R)) with multiplicity rrij. Then, by 
Lemma [TJ it follows, with m = max(m,;, rrij), that 

log \zi — Zj\ > -n 2 /m - (n 2 /m + Zn 2 j {2m i m i )) log(n 2 + 1) - 2n 2 / '{miVrij) log \R\ 

> — 16n 2 logn/m — n 2 / (miTrij) log |i?|oo = — 0{n r/m), 

(4.7) 

where we used that log \R\00 = C(log n ■ n(logn + r)). By the construction of /, 
Zi is the closest root to z. Hence, Lemma [2] yields 

log \R(z)\ > rrn log \z - Zi \ - 2{n 2 /m t - 1) log \R\ X - 2n 2 {l + 3/(2™,)) logn 

> m t \og{\ Zl - z,|/(1024n 6 )) - (2n 2 /m t ) • log^U - 5n 2 logn. 

(4.8) 

By combining (|4.7p and (|4.8[) . we obtain 

logLS(a) > log|.fl(z)| - 2n 2 = -0(m 4 ■ n 3 T/m l +n 3 r/m t +n 2 ) = -d(n 3 r). 

(4.9) 

In completely analogous manner, we show that \ogLB{(3) > — 0(n 3 r). 

Estimating the upper bounds. We now investigate in estimating the upper bounds 
UB(a,P,vS v) ) and UB(a,p,v^) for \u^\ and |«(*>| on A(a,/3). In order to do 
so, we apply Hadamard's inequality to the matrices and V^ y \ see Section[3] 
By analogy, these estimates then also extend to the upper bounds UB(a, (3, u^ x ') 
and UB(a^,v^) for |u<»| and \v^\ on A{a,0). 

In the realization, to compute the upper bounds on the cofactors, we first 
use interval arithmetic on a box in C 2 to estimate the absolute values of respec- 
tive matrix entries fZy and on A(a,0), and then apply Hadamard's bound. 
In the complexity analysis, we follow a slightly different but even simpler ap- 
proach: From the construction of A(a,f3), the disc A(a) has radius less than 



a (a, i?( y ))/4, and the disc A(f3) has radius less than cr(/3, R^)/A; see Section[3] 
for details. Hence, it is clear that both radii are bounded by B = 2 ( n ( T + lo s™), 
where B is the upper bound from (|4.4[) for the modulus of the roots of 
and R( y \ It follows that, for each point (x,y) G A(a,j3), we have \x\, \y\ < 2B. 
Recall that, the entries of are the coefficients f^ v \x) and g\ v \x) which are 
polynomials of degree n x < n and with integer coefficients of bitsize r or less. 
Thus, for (x, y) £ A(a, 0), we can bound them as follows: 

\f?\x)\ < £ l^l W < (n x + 1) • 1/^1(25)^ < 2" +1 n • rB\ 

3=0 

and a corresponding bound holds for Ifl^^x)! as well. Hence, it follows that 
log \U^\x, y)\ = 0(n\ogB + r + n) = C(n 2 (logn + r)). Now, by Hadamard's 

inequality \u^\ = \det(U^)\ < lla. where \\U^ v) \\ 2 is the 2-norm of 

the z-th row vector of U^ y \ Since the absolute value of each entry of EA») is 
bounded by 2°(' l2 ( lo s n+r )), we have log||C^ (j/) || 2 = 6{n 2 r). Taking the product 
of all bounds on the 2-norms yields log \vS v '\ = C(n 3 r) and, thus, 

UB(a, [3, u (y) ) = 6{n 3 T)). (4.10) 

This bound is also valid for UB{a, (3, v^), UB(a, &M x) ) and UB(a, /3, v^). 

The inclusion test Let us first assume that a candidate box B := B(a, (3) = 
1(a) x I(f3) C M 2 contains a real solution of ([l.ip . Using our bounds from (|4.9p 
and (|4.10p . the inequalities in (|3.5[) and (|3 .6[) rewrite as: 

\f(x ,y )\ + \g(x ,y )\ < 6 = 2~ & ^\ (4.11) 

where 5 is a certain threshold |/(xo, yo)\ + \g(xo, yo) \ has to undercut. In the case 
where B does not contain a real solution of (|1.1[) , Theorem [T] and (|4.11l) ensure 
that, for all (xo,yo) G B(a,f3), we have: 

\f(x a ,y )\ + \g(x 0l y )\>2- d ^l (4.12) 

Now, in order to certify or discard (a, (3) as a solution of f = g = 0, we set 
p := [~— log s] that is directly related to the size s of B and evaluate f(a, (3) and 
g(a,/3) by means of interval arithmetic with precision p; see Q3J Section 4] and 
Appendix IB1 for the definition of polynomial interval evaluation with precision p. 
Obviously, the so obtained intervals 9S(/(a, 13), p) and ( B(g(a, (3), p) then contain 
f(B) and g(B) because each (x , yo) G 5 approximates (a, (3) to an error 2~ p or 
less. Now, if one of the intervals 53(/(a, /3), p) or *B(.g(a, (3), p) does not contain 
zero, then (a, (3) cannot be solution of (|1.1[) . If, for all (ei,£2) G Q3(/(a, /3), p) x 
Q%(a, £),/>), both inequalities UB{a, (3, u^))|ei| + J7B(a, /?, i/ y) )|e 2 | < £-B(a) 
and [7B(o!,^,uW)|ei| + t/B(a, /3, v^>)|e 2 | < LB(/3) are fulfilled, then (a,/3) 
must be a solution. In the case where we cannot discard or certify (a,/3) as 



a solution, we refine B by refining the corresponding isolating intervals 1(a) 
and I (/3) and proceed. According to (j4.11[) and (|4.12p . we eventually succeed if 
Q3(/(a, (3), p) and *B(g(a, (3), p) have width less than some threshold larger than 
2-6 (n r) j n |_ nc f }i owm g consideration, we will bound the cost for the evalua- 
tion of <8(/(a, (3), p) and *B(/(a, /?), p) and for the refinement of B. 

Complexity of polynomial evaluation. In the previous section, we have shown 
that the candidate boxes B = B(a,f3) have to be refined to a size s such 
that s S(f(a, f3), p) and < 3(g(a, (3), p) have width less than a certain thresh- 
old larger than 2~ (™ r >, where we consider interval arithmetic with precision 
p = \— logs] . We first estimate the absolute error induced by the interval arith- 
metic. Then, we estimate the overall cost for the polynomial evaluations at all 
candidate solutions. 

Theorem 4. The cost for evaluating 23(/(a, (3), p) and *B(<?(a, (3), p) for all can- 
didate solutions (a, (3) € C using interval arithmetic with precision p € N is 
bounded by 0(n 5 (p + n 2 r)) bit operations. Each of the so obtained intervals has 
width less than 

2 -P 2 t B 0(«) _ 2 -p+6(n 2 r)^ 

Proof. Remark that evaluating a polynomial F of degree d at a point xq using 
fixed point arithmetic with precision p demands for 0(d(p + log(max{|a;o|, 1}))) 
bit operations. Since each root a of the resultant polynomial has modulus 
less than B = 2°( nT \ we can compute the coefficients fi(a) of f(a, y) (by means 
of interval arithmetic with precision p) using 0(n 2 (p + n\ogB)) bit operations, 
while the resulting intervals contain values of modulus less than 2 T ' B°( n \ Thus, 
the cost for evaluating all fiber polynomials f(a,y) is bounded by 0(n A (p + 
n\ogB)) because has at most n 2 many roots. By [14j Lem. 3], the absolute 
error for each coefficient fi(a) is given by 

\fi(a) - <B(fi(a),p)\ < 2-P +1 (n + l) 2 2 T B n = 2-n T B° {n \ (4.13) 

In a second step, we compute < S(f(a, /3), p) for a fixed a. Again, for each (3, 
we have to perform n multiplications and additions involving numbers of bit 
length 0(p + t + nlogB), hence, the cost for these evaluations is bounded by 
0(n(p+r+nlogB)) bit operations. Summing up over all (3 above a fixed a yields 
(D(n 3 (p+T+nlogB)) bit operations and, thus, 0(n 5 (p+T+nlogB)) = (D(n 5 (p+ 
n 2 r)) constitutes an upper bound on the number of bit operations needed for 
all evaluations. The absolute error for each f(a,(3) is bounded by 2~ p 2 r B°^ 
because each coefficient fi(a) is approximated to an error 

2 -p 2 r B 0(n) and the 

absolute value of each f3 is bounded by B. The corresponding bounds for the 
evaluation of *B(g(a, /3), p) follow in completely analogous manner. 

Now, according to the above Theorem and our considerations in the previous 
section, we must succeed for a precision p > p = <D(n 3 r). Namely, since the 
width of the intervals < B(f(a, f3), p) and < 3(g(a, f3),p) is bounded by 2-p +0 ^ t \ 
it suffices to consider a precision p = 0(n 3 r) to guarantee that the so obtained 



intervals have width less than a threshold that is lower bounded by 2 °( n T K 
Hence, the cost for all polynomial interval evaluations is bounded by 

6{n 5 (n 3 T + n 2 r)) = 6(n 8 r). 

Complexity of the root refinement. We have already shown in the previous section 
that it suffices to refine each of the candidate boxes B(a,0) = 1(a) x I((3) 
to a size 2~® < - n T \ As mentioned in Section [3j we use the Quadratic interval 
refinement (QIR) method to refine the isolating intervals 1(a) and 1(0) for the 
real roots a, (3 of R^ and R^ x \ respectively. Recent work [Tl] improves on the 
asymptotic complexity for this root refinement by considering approximations of 
polynomial coefficients to a certain precision. It is shown that, for a square- free 
polynomial F with magnitude (d, /j.) and root bound r, the algorithm refines 
isolating intervals for all real roots of F to a width 2~ L or less using 

6(d(d\o g r(F) + S(F))(fi + d\o g r(F) + 17(F)) + d 2 L) (4.14) 

bit operations, see [HI Thm. 22]. We need to refine the real roots of the square- 
free part R* of R = to a width bounded by 2~° ( " 3t) . Using the above 
result, we can do so with 

6(n 2 ■ n'V ■ (n(n + r) + n 3 r) + n 4 • n 3 r) = 6(n 8 T 2 ) (4.15) 

bit operations, where we used that S(R*) = (D(n 3 r) (see (|4.3p ) and all roots 
of RS V > have modulus less than B = 2 c, (™( log "+ r )) . Again, the corresponding 
bounds for refining the roots of R^ follow in complete analogous manner. 

5 Conclusions 

We have derived a bound on the bit complexity of the algorithm Bisolve to 
isolate the real solutions of a bivariate polynomial system. To the best of our 
knowledge, the derived complexity bound improves considerably upon the best 
known complexity bounds to date. Still, we suppose that the derived bounds 
are not tight and there is room for improvement. This mainly stems from the 
fact that the number of solutions of a bivariate system is in most cases largely 
overestimated while, at the same time, this number is intimately bound up with 
the root separation and coefficient bitlength of the resultant. We find it, there- 
fore, reasonable to invest some research in the theory of sparse resultants where 
the geometry of exponent vectors of input polynomials is analyzed to obtain 
sharp estimates. A bottleneck in the current analysis is the isolation of the re- 
sultant polynomial which is assumed to be a general polynomial of magnitude 
(n 2 ,n(r + logn)). Hence, the question arises whether isolating the roots of a 
resultant is possibly easier than of a general polynomial? Finally, it might be 
possible to improve the complexity of the validation step by considering asymp- 
totically fast algorithms for multipoint polynomial evaluation [15| . 
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A Proof of Theorem [H] 

We begin with the following auxiliary lemma. 

Lemma 3 (Davenport-Mahler). Let F(x) be a polynomial with integer coef- 
ficients of magnitude (d, fi) , and V be the set of distinct complex roots Z\ , . . . , Z r 
of F (r < d). If a directed graph Q = (V,E) satisfies conditions: (a) for every 
edge (a,/3) € E, \a\ < |/3|; (b) Q is acyclic; and (c) the in-degree of any node is 
at most 1, then 



For a proof of the above lemma, we refer to [9l Corollary 3.11] . It exploits the fact 
that we can always factorize F over Z and extract its square-free part F* € 1\x\. 
Then, the generalized Davenport-Mahler bound applied to F* automatically ex- 
tends to F because for F*, as the divisor of F, it holds that M(F*) < M(F)$ 
Now we are ready to show the main result. 

Proof of Theorem^ We begin with clustering the roots of F into subsets consist- 
ing of nearby roots. Without loss of generality, we assume that cr{z\, F) < ■ ■ ■ < 
a(z r , F). For ft £ N, we denote i(h) the maximal index i with a(zi, F) < 2~ h and 
U = U(h) := {z\, . . . , the corresponding set of roots Zi with a(zi,F) < 

2~ h . If h < log(l/cr(F)), then U contains at least two roots. Our goal is to 
partition U into disjoint subsets U\ 1 . . . ,Ui that contain closely located roots. It 
can be shown that for h < log(l/cr(F)), there exists such a partition of U that 
\Ui\ > 2 for alH = 1, . . . , I and \z-z'\ < d2~ h for all z, z' G U u see [H Lem. 19]. 

We consider a directed graph Qi on each Ui connecting consecutive roots of 
Ui in ascending order of their absolute values. Let Q := (U,E) be the union of 
all Qi. Remark that, Q is a directed graph on U, and the conditions of Lemma[3] 
are satisfied. Moreover, since each set Ui contains at least two roots, we must 
have #E > i(h)/2. Additionally, it holds that #E < i(h), where in corner case 
all roots of F belong to a single partition U\. Hence, by Lemma [3] 






and therefore, 



2d(/x + log(d+l)) 2d( A t + log(rf+l)) 
log 3 + log d + h h 



3 Here, M(F) denotes the Mahler measure for F. 



By combining the above inequality with h < log(l/cr(F)), and i(h) > 2 (because 
\Ui\ > 2), we conclude that: log(l/er(F)) < d(fj, + log(d + 1)) + 1. Otherwise, 
there would exist an h with i(h) < 2 which is not possible. For the bound on 
S(F), it suffices to consider only the roots Zi,...,Zk with separation < 1/2 
since all other roots contribute with at most d to the sum U(F). Let h max = 
\d(n + log(d + 1))] , then: 

S(F) =^2log(l/ a( Zi ,F)) < £ i(h) < dO*+Iog(d+l)) £ - = 0(df,log(dfi)). 

i=l h=l h=l 



B Approximate polynomial evaluation 

In this section we define some basic operations on approximate numbers and 
establish the error bounds for polynomial evaluation, we refer to [2J Section 4] 
for comprehensive discussion. For p G N and x £ K, we define: 

down(a;,p) = {k-2- p e R\k = [x-2 p \}, and up(x,p) = {k-2~ p £ R\k = \x-2 p ]}. 

Meaning that, x is included in an interval: Q5(a;,p) := [down(a;,p), wp(x,p)]. 
In what follows, we will omit p and write up(x) or Q3(:r) to simplify the nota- 
tion. Arithmetic operations on approximate numbers obey the rules of classical 
interval arithmetic. For two numbers x, y € K, we define: 

*B(x) + [down(x) + down(y), up(x) + up(y)], 

*B(x) — Q3(?/) :~ [down(x) — up(y),up(x) — down(y)], 



min {Hi^Hjiy)}, max {Hi(x)Hj(y)} 

ij={l,2} ij = {1,2} 

with H\{x) = down(x), and H2(x) = wp(x). 



Using the above rules, for a polynomial F and z€l, < B(F(z), p) can be defined 
by expanding the polynomial according to Horner scheme: 

«8(F(z)) = Q3(F ) + Q5(z) • (Q3(iq) + • (®(F 2 ) + ...)). 

The next lemma bounds the error of polynomial evaluation with precision p: 

Lemma 4. Let F be a polynomial with magnitude (d,p), c € R with \c\ < 2 V , 
and p £ N. Then, 

\F(c) - H(F(c), p)\ < 2- p+1 {d + i)V +cto , 

where H = {down, up}. In particular, <8(F{c),p) has width 2- p+2 (d+l) 2 2 p+dv 
or less. For a proof, see \14\ Lem. 3]. 

In essence, the lemma asserts that the absolute error result from approximate 
polynomial evaluation is linear in 2 _p , which has important consequences for us 
in Section |47~ 



